1 



On 



Analysis of Water Vapor spatio-temporal structure over the 
Madrid Area using GPS data 

G. Ruffini 1 , A. Rius 1 , L. Cucurull 1 , A. Flores 1 



Short title: ANALYSIS OF WV SPATIO-TEMPORAL STRUCTURE 



2 



Abstract. We have analyzed Zenith Wet Delay (ZWD) time series from an 
experiment over the Madrid (Spain) area obtained from 5 GPS receivers using two 
different techniques. In the first case a delay correlation analysis of the ZWD time-series 
has been carried out. We show that for this small network (with a spatial scale of 
less than 100 km) the correlation between the time series is very strong, and that 
using windowing techniques a reliable correlation delay time series can be produced for 
each pair of sites (10 such pairs are available). We use this delay time series together 
with a frozen flow model to estimate the velocity of a passing front, and compare the 
results to meteorological data and Numerical Weather Prediction output, showing good 
agreement. In the second approach, the data is analyzed using Empirical Orthogonal 
Functions. We demonstrate that the temporally demeaned and normalized analysis 
yields information about the passing of fronts, while the spatially demeaned data 
yields orographic information. A common second mode highlights the underlying wave 
behavior. 
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Introduction 

If an important goal for the GPS research community has been to test the limits of 
the geophysical measurement techniques derived GPS technology, a now pressing task 
is to use the newly available data for meteorological studies. We will focus here on 
trying to extract relevant information from the new type of data generated by the GPS 
measuring technique. In a previous publication we discussed the analysis of Zenith Wet 
Delay (ZWD) and gradients measured with GPS and Water Vapor Radiometers (WVR) 
\Ruffini et ai. 19~99\ . We will now analyze the spatio-temporal structure of the obtained 



GPS Zenith Wet Delay Time series. Similar studies have been carried out by [ Davit 



et ai., 1998j\ (although the scale of the network involved, the Swedish permanent GPS 



network, is significantly different, where it was already pointed out that GPS water 
vapor estimates can be very useful for studying the spatial progress of air masses. 

The refractivity of the neutral atmosphere at radio frequencies is given approximately 
by N « 77.6P/T + 3.73 x 10 5 P W /T 2 = N dr y + N wet , where P is the total pressure, 
P w is the water vapor partial pressure (both in mb), and T is the temperature (in 
K). The equivalent excess path length corresponding to a ray crossing the neutral 
atmosphere is given by \\Bevis et ai. 1991^ AL = 10~ 6 / N dl + S — Q. Here Q is the 



straight-line distance between satellite and receiver, and S is the geometric path length 
along the ray. In order to estimate Zenith Total Delay (ZTD) mapping functions are 
used, because GPS measurements are not, in general, in the zenith direction. Mapping 
functions model the dependence of the tropospheric delay on satellite elevation, making 
some assumptions about the tropospheric gas distribution. In the past, azimuthally 
symmetric models were employed, and the elevation dependent slant delay approximated 
by AL(e) « ALf y m dry (e) + AL™ et 

f^weti^): where vn we i (e) and m dry (e) are elevation 
mapping functions. Tropospheric delay gradient estimation is now possible and routinely 
carried out \\Bar- Sever et ai. 199$ \MacMiiian 1995 \. 



Two types of tropospheric gradients were considered in \Ei6segui et ai., 
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1998b; \Ruffini et at. 199% . Let the atmospheric refractivity be given by N(p,z) 



N (z) + V p N(z, p)\p =Q ■ P, where z is the height coordinate and p the horizontal 
displacement vector, and define V p N(z, p)\^ =0 = g(z), the horizontal refractivity 
gradient. On one hand, there is a gradient associated with azimuthal dependencies 
of delay observations at a GPS receiver. This GPS tropospheric slant-delay gradient 
(sd- gradient), G = (Gn,Ge), is defined by the non-azimuthally symmetric delay part in 
the GPS signal by AD(e,p) = m&(e) cot(e) G ■ p, where p = (cos 0, sin 0) = p/||p|| is the 
azimuth unit vector \ \Bar- Sever et al. 1998\ , \MacMillan 1993\ . On the other hand, Zq, 



the horizontal gradient estimated with the zenith delays from a network (zd- gradient) 
is more closely related to g. As discussed in \ El6segui et al., 19981 ], one can show 
G = 10- e I^dzzg(z), Z = lO- e f™dzg(z). 

Here, we analyze the data from a GPS campaign carried out in the Madrid area 
during December 1996 in the light of a frozen flow model. It is known from analysis 
of meteorological data that two humid, cold fronts crossed the network during the 
campaign \\Cucurull et al. 1998^ . See also Figure 2 below. 



The basic idea in the frozen flow model is that the wet refractivity field propagates 
like a wave in the presence of a passing front: 

N((,p,t) = N Q (Z-p-ut).e-V h , (1) 

where the vectors are (2D) surface vectors. Here ( is a generalized vertical coordinate 
that depends on the distance from the geoid (z) but which may also depend on the 
orography, C, = z + i](p), and h is a vertical scale. The velocity associated with this wave 
is given by v = ku/n. Notice that temporal and spatial gradients are closely related in 
this model: v ■ V p N = V t N - NV p n/h or, equivalently, v ■ VW = f t W - WV p r]/h, if 
the ZWD (denoted by W) is measured at constant z — not the present case. 

Notice also that the refractivity gradient is closely related to k if we assume a 
geopotentially stratified atmosphere, i.e., if z = (, since by equation [l], we have would 
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have g = V p N(z, p)\$ =0 = HN^—ut^e^^ — N(—ut))V p r)/h, and the last term would 
drop out. Unfortunately, this approximation was not correct for our network, where 
orography plated an important role. 



GPS Campaign 

In the GPS campaign (December 1996), we deployed 5 Trimble geodetic GPS 
receiver systems (called ROBL, ESCO, VILA, IGNE and VALD) near the Madrid area 
Spain, on December 2-15, 1996, with inter-site separations from 5 to 50 km \ \hJlosegu 



\et at., 1998a\\ (see Table ). GPS observations consisted of data streams of undifferenced 
dual-frequency carrier-phase and pseudo-range measurements obtained every 30 seconds. 
The GIPSY/OASIS-II (v.4) software package || Webb et at. 199% (Gipsy) has been 
used with a Point Positioning strategy to estimate ZTD at the five GPS sites with a 
precision of 5 mm. Estimates of the satellite clock corrections and orbits were provided 
by the IGS and JPL, as well as consistent earth- rotation parameters. Gipsy uses a 
kalman filtering technique to model time-dependent observables, such as the ZTD. The 
tropospheric delay was modeled as a random walk, a 2 = d 2 ■ t, with a drift rate of d = 
0.25 cm/v^h. The drift rate for the gradient parameters was 0.03 cm/y/h. We used a 
cut-off elevation angle of 7° (see \\Ruffini et at. 1999} for more data processing details). 
The dry part of the delay can be estimated well (to less than 0.35 cm in delay) if surface 
pressure is known to within 1.5 mb. We have used pressure estimates produced by 
HIRLAM together with ground measurements, since barometric measurements were not 
available at all sites. A conservative estimate of the accuracy of the pressure data is 



better than a 1.5 mb \Cucurutt et at. 199 1 
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Correlation analysis 



Given the frozen flow model in equation [I], the time series of ZWD at the different 



sites can be approximated by 



Wi(t)=f{t + n)e-«, 



(2) 



where Wi(t) represents the ZWD at the z'-th site at time t, Tj represents a delay relative 
to some chosen site, and is the scaling factor that should depend partly on the 
height of the site (without loss of generality, we refer this height to VILA's). Thus, the 
delay Ty between two sites will be given by = ■ k = Xij ■ v/v, where the inverse 
velocity vector k is given by v/v, the quotient of the velocity unit vector and its norm, 
and Tij = Ti — Tj. The delay is calculated by finding where the expression 
F[rij] = (Wi(t)Wj(t + Tij))A-i attains its maximum. The t-superscript indicates that the 
time series' have been temporally demeaned and normalized to unit standard deviation, 



W l = {W — (W) t )/a t , where a t = y ((W — (W) t ) 2 )t- Cross-correlations are found using 

a time window given by A (twelve hours were used here). In Figure 3 we see an example 
of the analysis, where the signal is clearly seen. Note the presence of intervals where the 
correlation is high, and the delay is zero, which will be discussed below. 

The inverse velocity vector k has then been estimated, once for each time, by 
minimizing 



The results are plotted in Figure 4, and will be discussed below. In particular, we 
will try to interpret the existence of intervals where the estimated velocity is infinite, 
associated with zero delays. 





(3) 



EOF analysis 

Empirical Orthogonal Function (EOF) decomposition is a standard tool in 
multi-variate data analysis. For the task at hand, it will be convenient to think of the 
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ZWD time series at the different sites as a time series of images, I(t), representing at 
each time the WV content over a 2D network. Given these time series, the goal of EOF 
analysis is to decompose them in orthogonal modes, that is 

I(t) = j:\ 1 IyuJ(t). (4) 



We have used Singular Value Decomposition, a very useful tool in this context \\Keine7 
et ai, 199% 



Carrying out EOF analysis with the original time series, however, yields a very 
strong first mode (Xj = 338.9,8.2,5.6,3.0,4.2) which essentially represents the mean 
temporal behavior and the exponential vertical behavior. The two effects, however, are 
hard to separate and interpret, as spatial and temporal effects tend to get mixed up in 
the modes. For this reason it is useful to demean and normalize the data first, as we 
discuss below. 

It is important to emphasize that the most difficult part in EOF analysis is not the 
numerical computation of the modes (that is actually fairly simple), but to interpret the 
results. For this task we will make extensive use of the flow model. 

Temporal demeaning 

Let us first study the case in which we temporally demean and normalize the time 
series, since this parallels the correlation approach in the previous section. Recall that 
we model by equation |2|. It is the straightforward to show that 

w u t) _ fit + Tj) - (/(*)>* M m - wh | f^ T \ (5) 

St St St 

where sf = \(f(t) — (f(t)) t ) 2 } — assuming (f(t)) t ~ 0. Intuitively, EOF analysis in 
this case should yield information analogous to that of correlation analysis. We should 
obtain a spatially homogeneous first mode, with a strong temporal variation, and a 
second mode with a spatial structure related to the passing front (again, through the 
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delays involved). This spatial structure is directly related to the zd-gradients defined 
above, except for orographic corrections. 

Another consequence of equation |5] is that W-(t) = g(t + Tj). This means, for 
example, in the case of two GPS sites that the observation matrix will have two columns 
as follows, A = [g(t) g(t + r)]. Carrying out the Singular Value Decomposition yields 
A = UWV T with eigenvalues w± = Jl± (g(t)g(t + r)) t , corresponding V eigenvectors 
v± = (1, ±l)/\/2, and U eigenvectors u± = (g(t) ± g(t + t))/w±. Thus, we see that the 
second mode carries information about the spatial relationships between the delays. The 
second spatial eigenvector represents a spatial derivative, while the temporal eigenvector 
is akin to a time derivative. 

The case with three functions is more difficult to analyze. Let gi(t) = g(t + n), and 
9ij = (9i{t)9j(t))t- Let also s = ^ g ^ +9 ^ + l2i. ; anc i / — g\292z9z\- Then, the eigenvalues 
are given by A = 1 + 2s (cos^ 1 {l/s 3 ) + n2ir/3) , with n — 0, 1, 2. For example, if 
9i (t) = e-(*-^) 2 / 2A 7AV4^, we find (g t (t)g 3 (t)) t = e~ T ^ 2 . 

See Figure 5 for a representation of the first three modes using GPS data for one of 
the periods which involve the passage of a front (the eigenvalues for the decomposition 
were A = 10.0,3.3,2.3, 1.4, 1.2). Our simulations of passing fronts yield precisely this 
structure. We have generated time series simply by taking one of the real ones and 
shifting it in time in a manner conforming to that of a passing front. An example of the 
resulting EOF analysis is shown in Figure 6 — a north directed front was simulated with 
a speed of 60 km/h. The eigenvalues in this case are A = 12.37,1.6,0.6,0.4,0.0. The 
fact that the last eigenvalue is zero is due to the need to specify only 4 numbers in the 
simulation (the Tj's). 

Spatial demeaning 

In this case, spatial demeaning and normalization to unit variance are carried out 
prior the EOF decomposition, W x = (W - (W) x )/a x , with a x = yf ((W - (W) x ) 2 ) x . 
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The x subscript means that statistics are to be computed using the site index i. From 
equation 0, we find, 

wx „ e~<* - (e<% d(\nf(t)) r t e^ - {r^% 

a dt s ' 1 ' 



where s x = y {(e~^/ h — (e~&/ h ) x ) 2 ) x . Hence, up to normalization, we may expect modes 
such as voi ~ — (e - "), with very little temporal variation, and vu ~ Xje~^ — (rje - ^), 
with a temporal variation associated again to the time derivative of f(t). The first 
mode should be a reflection of the orography of the network, and the next modes 
should encode the delay structure associated with the passing front. See Figure 7 for an 
illustration of the first 3 modes. The resulting eigenvalues were A = 0.6, 0.3, 0.2, 0.1, 0.0. 
Note here the disappearance of the last mode. 

Spectral analysis 

The spatially interpolated time-series can also be analyzed spectrally. This analysis, 
however, is not simple to carry out or interpret, as we will see. 

Spectral anaysis of the second EOF mode between the times of 13.6 and 14.2 
days reveals peaks at the harmonic spatial frequency with several peaks in temporal 
frequencies, including one at 0.7 per hour. This leads to a north-east velocities of around 
40 km/h. The measured average surface wind speed was actually 25 km/h, and wind 
direction was 42 degrees. 

What should we expect from this analysis? A simple model for a traveling wave is 
given by a gaussian wave-packet, Q = e -( kx ~ wt ) 2 . The Fourier transform is given by 



F(k, u) = ^e- k2 /* k ' 2 5{uj'/k' - u/k). (7) 
kk' 



For a more general wave in a non-dispersive medium, Q = Q(k ■ x — wt), we obtain a 
similar result. The salient feature is a diagonal spectrum along the constant velocity 
line k/u = k'/u'. In 2+1 dimensions the result is 

F{k, u) = 5 2 {k/u - k'/u') f e iu)u Q{u)du. (8 
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As may be expected, perfect simulations with the required spatial density illustrate 
this behavior, yielding a strong diagonal feature. Using the simulated data mentioned 
above, however, we find a north-direction-time spectrum with an axis of symmetry 
defined by a constant velocity vector, but not a diagonal spectrum. This is because the 
simulation, although perfect at the station sites, loses coherence when interpolated in a 
uniform grid, and the diagonal feature seems to be very unstable. The east direction 
does conform to the k e = equation. The resulting spectra for the north simulation 
and the real time series are very similar, however, including the k e = equation. We 
conclude that spectral analysis of the raw time series are of limited use. 

Conclusions 

We have shown here that cross- correlation and EOF analysis can be very useful 
tools for the detection of passing wet/dry fronts in small to medium size networks. 

A striking feature of the EOF analysis is the similarity in the secondary modes 
in the spatial and temporal demeaning approaches. As has been discussed before (see 
Keiner et al., 199% and references therein), this is a feature of modes whose temporal 



oscillations create spatial gradients, as is the case in the frozen flow model — a simple 
wave model. 

We can also compare our results with those in \\Rujjini et al. 1999\ . We can see 



there the passage of the front detected here in the form of the obtained sd-gradients 
(see Figure 8), although it is hard at this point to make very quantitative statements. 
The signal seems to appear in the two analysis, however, as we can see a south pointing 
gradient. This in, in effect, a comparison of zd-gradients with sd-gradients. If the 
exponential law were exact we would see of course a match between the two. This is not 
the case, however, because orography plays an important role. 

We hypothesize that infinite propagation speeds associated the zero delay 
correlations (see Figure 1) are related to rain events. For instance, we can imagine that 
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at some time it is raining at all sites simultaneously: a drop in WV will be recorded at 
all sites, yielding high correlation with zero delays. Any phenomenon that can change 
the WV content estimates in the network at the same time will produce this effect, 
however, so the conclusion is not warranted. We have plotted rain rate measurements 
in the Figure as well, for comparison. The peak of rain rate does seem to be associated 
with the zero delay period and a ZWD drop during that time. 
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Table 1. Positions of the re- 
ceivers with respect to VILA (all 
in km except AXz, in meters. 



site 


AX N 


AX E 


AX Z 


ESCO 


15.964 


-16.834 


430 


IGNE 


0.279 


20.576 


119 


ROBL 


-1.558 


-25.256 


181 


VALD 


4.852 


-7.158 


197 


VILA 


0.000 


0.000 
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Figure 1. The ZWD time series for the GPS and WVR during the campaign. The 
periods during which no WVR data is available are associated with rain events, as the 
WVR cannot work in wet conditions. 
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Figure 2. A Meteosat IR photo for December 12 1996 noon UTC (MET5 12 DEC 1996 

1200 IR1 D2). 
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Figure 3. The delay structure between ESCO and ROBL (correlation is shown dashed 
and normalized to 100 for graphing purposes). The fact that the delays are negative 
mean, in the convention used, that ROBL detected the changes before ESCO. In the 
bottom panel ROBL is shown dashed. 
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Figure 4. Results of the correlation analysis, and their comparison with meteorological 
data. In the top panel we see the ZWD at VILA, and the average correlation function 
of the analysis (normalized to 10). In the second panel we see the station mean delay 
from the correlation analysis and the rain rate. In the next panel the estimated k is 
plotted. In the panel below the ground mean speed (in blue) is plotted vs. the estimated 
wave speed. Finally, the estimated wave direction is plotted against the measured wind 
direction (in blue) in the bottom panel. 
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Figure 5. The first three modes after temporal demeaning. A = 10.0, 3.3, 2.3, 1.2, 1.4. 




Figure 6. The first three modes after temporal demeaning with a simulated 60 km/h 
north going front. A = 12.37, 1.6, 0.6, 0.4, 0.0. 
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Figure 7. The first three modes after spatial demeaning. A = 0.6,0.3,0.2,0.1,0.0. 
The first mode is closely associated to the orography. 
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Figure 8. Gradients obtained with GPS: VILA (solid), ESCO (dots), ROBL (dashes). 



